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Abstract 

The potentially significant role of the surface of an elastic body in the overall response of the 
continuum can be described using the mature theory of surface elasticity. The objective of this 
contribution is to detail the finite element approximation of the underlying governing equations 
(both in the volume and on its surface) and their solution using the open-source finite element 
library deal. II. The fully-nonlinear (geometric and material) setting is considered. The non¬ 
linear problem is solved using a Newton-Raphson procedure wherein the tangent contributions 
from the volume and surface are computed exactly. The finite element formulation is imple¬ 
mented within the total Lagrangian framework and a Bubnov-Galerkin spatial discretization of 
the volume and the surface employed. The surface is assumed material. A map between the 
degrees of freedom on the surface and on the boundary of the volume is used to allocate the 
contribution from the surface to the global system matrix and residual vector. The deal. II 
library greatly facilitates the computation of the various surface operators, allowing the numer¬ 
ical implementation to closely match the theory developed in a companion paper. Key features 
of the theory and the numerical implementation are elucidated using a series of benchmark 
example problems. The full, documented source code is provided. 

1 Introduction 

The surface elasticity theory of Gurtin and Murdoch (1975) has been widely used to account for 
the role that the surface of an elastic body can play in the overall response of the continuum. 
Integral to the theory is the derivation of a set of governing equations and constitutive relations 
that describe the behaviour of the surface of the bulk object. The role of surface elasticity and the 
size-dependence of the elastic response has received considerable attention recently (see e.g. Duan 
et ah, 2009; Weissmuller et ah, 2010). This resurgence of interest in the mechanics of surfaces can 
be largely attributed to the increasing number of applications involving nanoscale structures. In 
such applications scale effects are observed due to the significant surface-to-volume ratio. Classical 
continuum formulations are unable to account for these effects as they lack an inherent length scale. 

The theory of surface elasticity is well understood (see e.g. the review by Javili et ah, 2013a). 
The vast majority of numerical treatments of surface elasticity, however, have been limited to the 
infinitesimal deformation regime. Furthermore, to the best of the authors knowledge, no open-source 
or commercial finite element implementation of surface elasticity is available. 

The objective of this contribution is to detail a finite element implementation of surface elasticity 
at finite deformations using the open-source library deal. II (Bangerth et ah, 2007, 2013a,b). The 
implementation uses various deal. II routines developed to facilitate solving partial differential 
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equations on curved manifolds (see e.g. DeSimone et al., 2009; Heltai, 2008) embedded in a higher¬ 
dimensional space. The deal. II term for this set of routines is codimension one. The name follows 
from the mathematical definition of codimension: if IT is a subspace of a linear space or manifold 
T, then the codimension of IT G T is defined by codim IT = dim V — dim IT. The codimension 
one routines greatly simplify the construction of the approximations to the various mathematical 
operators on the surface, and the geometric description of the surface. This is critical, as the surface 
deformation gradient (used to parameterize the constitutive response) is rank deficient, and the 
surface profile can be complex when modelling realistic problems. The codimension one routines 
allow the surface of the body to be treated as an independent two-dimensional manifold embedded 
in three-dimensional space. The constitutive model for the surface implemented here allows the 
surface to behave in a solid- or fiuid-hke manner. For solid-like behaviour, the surface free energy 
resembles a classical neo-Hookean model. The neo-Hookean model is extended to include surface 
tension and thereby account for fluid-like effects. 

Certain aspects of the implementation presented here for the volume contributions are similar to 
those discussed in the online tutorial (step_44) on a three-field formulation for (near-incompressible) 
finite elasticity (see www.deaIii.org/deveIoper/doxygen/deaI.II/step_44.html). 

The finite element formulation is implemented within the total Lagrangian framework and a 
Bubnov-Galerkin spatial discretization of the volume and surface employed. The surface is assumed 
material (i.e. it acts like a membrane permanently attached to the underlying solid volume). The 
nonlinear problem is solved using a Newton-Raphson procedure wherein the tangent contributions 
from both the volume and surface are computed exactly. A curvilinear-coordinate-based finite 
element methodology for the problem of surface elasticity was developed in a companion paper 
(Javili et ah, 2013b). Extensive theoretical details and references are provided. For additional 
information on the admissible range for the surface material parameters, see Javili et al. (2012). In 
particular, the validity of negative surface parameters, which have been reported in the literature, 
is assessed. 

While the majority of the numerical examples presented in the literature are performed on 
academic problems, one intended and novel application of the code presented here is to describe 
realistic nanostructures with complex surface geometry. This requires an efficient and parallelized 
code. In order to realise this objective, various numerical operations in the volume and on the surface 
are parallelized (within a shared memory paradigm) using the Threading Building Blocks library 
(TBB) (TBB, 2013). Aspects of the extension of this framework to a distributed environment is 
discussed briefly. 

The structure of the paper is as follows. The key kinematic concepts are recalled in Section 2.1. 
Thereafter, the independent hyperelastic constitutive relations governing the response of both the 
volume and the surface are summarised. The strong form of the governing equations and their 
restatement in weak form as a Newton scheme are presented in Section 2.3. The linearized weak form 
provides the basis for the fully-discrete problem introduced in Section 3. Details of the numerical 
implementation within the deal. II library are given in Section 4. Four numerical examples are 
presented in Section 5. They serve to elucidate key features of the theory and aspects of the 
numerical implementation. The complete, documented source code, instructions, and all input 
decks required to reproduce the numerical examples are provided online at www. cerecam.uct. ac. 
za/code/surf ace_energy/doc/html/index.html. 

Notation and definitions 

Direct notation is adopted throughout. Occasional use is made of index notation, the summation 
convention for repeated indices being implied. When the repeated indices are lower-case italic letters, 
the summation is over the range {1,2,3}. If they are lower-case Greek letters the summation is over 
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the range {1,2}. The scalar product of two vectors a and b is denoted a b = [a]f[6]i. The scalar 
product of two second-order tensors A and B is denoted A : B = [A]ij[B]ij. The composition 
of two second-order tensors A and B^ denoted A B^ is a second-order tensor with components 
■ B]ij — mj- The tensor product of two vectors a and 6 is a second-order tensor 

D = a0b with [D]ij = The two non-standard tensor products of two second-order tensors 

A and B are the fourth-order tensors [A^ B]ijki = [A]ik[B]ji and [A0B]ijki = [A]ii[B]jk. 

An arbitrary quantity in the volume is denoted as {•} and analogously {•} denotes an arbitrary 
surface quantity. The surface quantity can be a vector, not necessarily tangent to the surface, or a 
tensor, not necessarily tangential or superficial to the surface. The (conventional) identity tensor in 
is denoted as i in the spatial configuration and I in the material configuration. In what follows 
the identity tensors i and I are understood as the conventional identity tensors in E^, i.e. their 
matrix representation would be a 3 x 3 matrix with 1 in the main diagonal entries and 0 elsewhere. 
Although these identity tensors are invariant and z = J, we use different letters to indicate explicitly 
which configuration they belong to. 


2 Summary of the problem of surface elasticity 

The objective of this section is to recall briefly the problem of surface elasticity at finite deformations. 
For additional details the reader is referred to Gurtin and Murdoch (1975); Javili et al. (2013a,b), 
among others. 

2.1 Kinematics 

Consider a continuum body that takes the material configuration Bq at time t = to and is mapped 
via the (volume) non-linear deformation map cp to the spatial configuration Bt at time t > 0, as 
shown in Fig. 1. Material points in the material and spatial configurations are denoted X and 
X, respectively. The associated linear deformation map, i.e. the (volume) deformation gradient, is 
denoted hy F := dxjdX = g-0G\ and maps material line elements dX G TSq (tangent to So) to 
spatial line elements dx ^ TBt via the relation dx = F -dX. The co- and contravariant basis vectors 
in the material and spatial configurations are denoted Gi and respectively. The displacement of 
a material point is denoted u = p>{X, t) — X. The volume deformation gradient F is rank-sufficient 
with inverse / := dXjdx = Gi g'^. The determinant of the deformation gradient and its inverse 
are given by J = det F > 0 and j = det / = 1/J > 0, respectively. 

Let So and St denote the surface of the continuum body in the material and spatial configurations, 
respectively. The outward unit normals to So and St are denoted N and n, respectively. Material 
particles on the surface are denoted X in the material configuration and are attached to the volume, 
i.e. X = X 1^230- Consequently, So = OBq. Furthermore, we assume that the surface is material 
in the sense that it is permanently attached to the substrate (i.e. the boundary of the volume). 
Therefore St = dBt and x = x\Qts^. This assumption implies that the motion of the surface tp is the 
restriction of the volume motion cp to the surface, i.e. ^ = (pldBo- 

Material line elements on the surface in the material and spatial configurations are denoted 
dX G TSo and dx G TSt, respectively, and are related as dS = F • dX, where F := dxjdX = 
g^ (g) G^ denotes the rank-deficient, and thus non-invertible, surface deformation gradient. The co- 
and contravariant surface basis vectors in the material and spatial configurations are denoted G^ 
and respectively. Various key concepts from differential geometry are given in Summary 1. The 
inverse deformation gradient / is defined by / := dX jdx = Ga C and is related to the surface 


3 



and 


deformation gradient as follows: 

f = I =: I - N 0N 


F'f = i=:i — n0n. 


The (surface) determinant of the surface deformation gradient and its inverse are denoted J = 
detF > 0 and j = det/ > 0 , respectively, where j J. 



Figure 1: The material and spatial configurations of a continuum body and its boundary and the 
associated deformation maps and deformation gradients. 


2.2 Constitutive relations 

A hyperelastic, neo-Hookean material model is assumed for the constitutive response of the volume 
(see Summary 2 ). The free energy T = T(JP) is parametrised by the Lame moduli A and ji. The 
constitutive relation on the surface is chosen to mimic that in the volume. In addition to the neo- 
Hookean type surface hyperelastic response, the surface free energy T = T(F) accounts for surface 
tension via the parameter 7 . The contribution of surface tension renders the surface free energy 
non-zero in the material configuration. This has implications for the resulting numerical scheme. 
For a study of the admissible range for the surface material parameters, see Javili et al. ( 2012 ). 

2.3 Governing equations 

The strong form of the equations governing the response of the volume and the surface (i.e. the 
quasi-static balances of linear momentum) are given by 

DivP 4 - 6 q = 0 in So , (1) 

DivP — P ■ N = 0 on So , ( 2 ) 

where 6 q is the (prescribed) body force per unit reference volume in the material configuration, 
and 60 is the (prescribed) force per unit reference area of the material configuration. Note that in 
the absence of a surface with an independent free energy (i.e. an energetic surface), ( 2 ) defines the 
standard Neumann boundary condition on the Piola traction, i.e. = P • N. 

The weak form is obtained by respectively testing (1) and (2) with arbitrary test functions 
Sep G Hq{Bo)^ and Sep G iLo(So)^, applying the divergence theorem in its extended form (see e.g. 
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Summary 1 The key differential geometry concepts of the surface (see Javili et ah, 2013b). The 
surface coordinates are denoted 


dr = dr(|) = dr(^^,(^^) 



^ dr 

9^ =91^92 , 93 


gradji} = 0 

d£,°‘ 


det{»} 


r = ^ with aG{l,2} 

, n = 

’ div{*} = = grad{i} : i 

d£,°‘ 

|[{»}-gi]x[{»}-g2]| 

\9i X 92\ 


9a=9a0 9'^, 9a0=9a-90, 9“ = 9/3 , 9°"^ =9°" '9^ , [fla/s] = [9"'^] \ 9 = |[?a/3]| 

g : surface permutation tensor, g = 9a/39“ 0 9^ 0 9^ = 9a ® 9/3 ® 93 > 9a0 =9°^^ 9 = \/d ea/3 


1 if ap is 12 

ea/3 = •{ -1 if a/3 is 21 , ?a /3 = | [ 9 a x 9 / 3 ]I = \/? 6 ^/ 3 , 9 “'^ = |[9“ X g^]| = [y^]-! ea /3 

0 otherwise 


uxv = [u^v]\^,U'V = [u^v]\i,i = 5^g^^g^=g^ 


gi^g^^g2^d‘^ = 'i'-^^T^ 


Javili et ah, 2013a), exploiting the orthogonality properties of the surface Piola-Kirchhoff stress 
measure (i.e. P • N = 0), and the kinematic constraint on the motion of the surface dip = 5(p\dBoi 
rendering the weak form as: 



where the dependence of the Piola-Kirchhoff stress measures P and P on the solution is given via 
the constitutive relations in Summary 2. The Neumann part of the surface is denoted Sq . Note 
that from the assumption of a material surface the final term of (3) is thus zero due 

to traction continuity (as indicated). It will prove convenient, when constructing the finite element 
approximation, to decompose (3) into volume and surface contributions. The volume and surface 
contributions will be recombined when the linear system is assembled. 

The decomposed weak form of the governing equations is solved using a Newton-Raphson strat¬ 
egy. The time domain [0,T] is decomposed into N uniform intervals of duration At := T/N = 
tn+i — tn^ where tn = nAt. The complete state of the system is assumed known at tn- The (not 
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necessarily converged) value of a variable evaluated at an iteration (i) during the n + 1 timestep is 
denoted The volume and surface residual contributions, denoted H and H, respec¬ 

tively, and the directional derivatives of the residual contributions in the direction of the solution 
increment in the volume Au = and on the surface Au are defined by 


and 



-J 

Grad^(^ 


Bo 



-J 

Grad^^ 


So 



■-j 

Grad^(/:? 


Bo 



■-J 

Grad^^ 


<5o 


dy 


p(i) dA 


- J 6<fi-bPdV, 

J3o 

- I Slp-P^dA, 


dp(i) 

dF 


Grad Am dV, 


dpd) -_ 

—^ : GradAndA. 

OF 


The derivatives of the Piola-Kirchhoff stress measures with respect to the deformation measures (i.e. 
the tangents) are given in Summary 2. Note, an exact computation of the tangent in the volume 
and on the surface is used. 

The resulting Newton scheme is thus given by 


pd + + Rii) + = 0 , 


+ Dag.R(*) = -Pd _ . 


(4) 


3 The fully-discrete problem 

The material volume and the surface are partitioned into sets of non-overlapping cells (elements) 
individually denoted f7e,o and 5 respectively (see Fig. 2). The discrete form of the governing 
equations is obtained by approximating the test functions and trial solutions in the linearized weak 
form (4) using vector-valued shape functions. The volume and surface shape functions, associated 
with an arbitrary degree of freedom in the volume I and on the surface /, are denoted and 
respectively. The associated value of the degree of freedom (i.e. the increment in the motion) in the 
volume and on the surface are denoted and respectively. They form the entries of the global 
solution vectors u and u. The trial solution in the volume and on the surface, and their gradients 


are approximated 

as follows^ 




u ^ 

J 

and 

Gradn 

a y]Grad#'^u^ 

J 

u ^ 

: y] 

and 

GradS 



I I 


The same approximations are used for the test functions (i.e. a Bubnov-Galerkin spatial discretiza¬ 
tion is employed). 

^Here and henceforth we adopt the abbreviated summation notation: where n^of is the number of 

degrees of freedom in the volume (or on the surface). 
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volume 


surface 



Figure 2: The triangulation of the material volume Th and a typical volume cell r^e,o- The surface 
triangulation Th is extracted from the volume triangulation. A typical cell on the surface is denoted 

^e,0 • 


The resulting linearised incremental problem to be solved is given in matrix form as 

AW + A«1 [Au] = - [rW + R(i) 


(5) 


□ □ 

where the global residual vector R and system matrix A are assembled from the element contributions 
as follows: 


where 


rw = 

^ I 

rw = Ay^R^ 


and 

and 


^ I J 

aw = AeEaF. 


/ J 


/ 


R^ := / Grad#^:Pdy- / 


/ 


rI:= j Grad$^:PdA- j 

^ e ,0 ^^,0 

Ay := j Grad#^ : A : Grad#-^ dF, 


Ay := j Grad$^ : A : Gr^^'^dA , 


and A denotes the (non-standard) assembly operator. The assembly operator is defined such that 

e 

contributions from degrees of freedom on the surface are mapped to the corresponding degree of 
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freedom in the volume (the surface is material and the triangulation of the surfaces matches the 
triangulation of the boundary of the volume). The number of surface degrees of freedom will always 
be less than those in the volume. The assembly operator is defined such that the size of the matrices 
and are the same. In other words A^*^ will contain empty rows and columns corresponding 
to degrees of freedom internal to the volume. 

4 Numerical implementation using the finite element library 
deal.II 

The open-source finite element library deal. II (Bangerth et ah, 2007, 2013b) is used to assemble 
and solve the discrete system of equations governing the problem of surface elasticity. The im¬ 
plementation here uses various deal. II routines developed for the solution of partial differential 
equations on curved manifolds (see e.g. DeSimone et ah, 2009; Heltai, 2008). 

The objective of this section is to review various key features of the numerical implementation, 
and to discuss some of the overarching design decisions. The complete, documented source-code and 
instructions can be found online at www.cerecain.uct.ac.za/code/surface_energy/doc/html/ 
index.html. Certain aspects of the implementation are similar to those discussed in the online 
tutorial (step_44) on (near-incompressible) finite elasticity without surface effects (see www.dealii. 
org/developer/doxygen/deal.II/step_44.html). 

4.1 Triangulations and degree of freedom handlers 

The triangulation (a class in deal. II) of the volume, denoted 7^, consists of the location of 
the vertices of the volume cells (elements) and the cell connectivity, as depicted in Fig. 2. A 
degree of freedom handler (another class in deal. II) combines the purely geometrical informa¬ 
tion of the triangulation with details of the finite element interpolation space. The surface de¬ 
gree of freedom handler is extracted directly from the volume mesh using the function 
GridTools: : extract_boundary_mesh. The surface degree of freedom handler is independent to 
that of the volume. A map Al, named surf ace_to_volume_dof _map in the code, is used to link the 
degrees of freedom on the surface to those in the volume. The surface triangulation is denoted 77- 

The contribution from the volume to the global system matrix is obtained by looping over all the 
volume cells f2e,0 5 assembling the matrix contribution from the cell Ag, and then adding this to the 
global system matrix. This is governed by the method SoIid< assembIe_system_voIume. The 

contribution from the surface is obtained in a similar way. A loop is performed over all surface cells 
Qe,o and the tangent contribution Ag calculated. The map A4 between the degrees of freedom on the 
surface and the corresponding ones in the volume is then used to add the surface cell contributions 
to the global system matrix (recall the surface is material). 

The Newton scheme (5) continues until the normalised measure (relative to the first iteration 

□ 

of the current time step) of the solution increment Au and the residual R decreases below a user- 
specified tolerance. 

It is important to note that alternatively a separate global surface system matrix could be 
constructed as the degrees of freedom on the surface are independent to those of the volume. The 
degrees of freedom on the surface could then be constrained to be the same as those in the volume 
and the resulting global block system solved. Using this approach, one could also then consider 
the case of non-material surfaces. It was decided that this approach would be more cumbersome 
as the global system matrix would have to be initialized using sparsity patterns from different 
triangulations. The approach adopted here can still be used to solve membrane problems (i.e. where 
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the volume is absent) by defining a volume with vanishing strength. More details are given in the 
liquid bridge numerical example discussed in Sect. 5.3. It should be emphasised, however, that the 
solution of membrane problems is not the focus of this contribution. 

4.2 The ContinuumPoint class 

The ContinuumPoint class represents a Lagrangian continuum point in the volume or on the surface. 
Its primary role is to calculate the constitutive response at a quadrature point given the kinematic 
state. For example, it determines the Piola-Kirchhoff stress and the tangent based on the defor¬ 
mation gradient, and its determinant and inverse. It is also convenient to store the kinetic and 
kinematic data in the class structure for use when constructing the residual and for post-processing 
the results. One key design decision was not to create separate classes for the material response and 
the quadrature point data as done in step_44 as this leads to data ownership conflicts. 

An instance of a ContinuumPoint is aware if it is in the volume or on the surface. Such knowledge 
is required occasionally (for example, when computing the tangent) but for the majority of the 
operations it is not required (for example, when computing the stress). 

4.3 Computation of kinematic quantities 

The form of the constitutive response at a quadrature point within the volume is near identical to 
that on the surface (see Summary 2). The primary difference is the structure of the deformation 
gradient (and its inverse) used to parametrise the free energy. The computation of the deformation 
gradient in the volume and on the surface are performed here in an identical manner (as depicted 
in Fig. 3) using the routines Solid<spacedim>: :update_volume_cp_incremental_one_cell and 
Solid<spacedim>::update_surface_cp_incremental_one_cell. 

The reference volume cell (i.e. the isoparametric domain) is denoted Cl\j. A typical volume cell in 
the material configuration is denoted ^4e,o- A spatial view of the volume cell, denoted is obtained 
by applying the nonlinear motion map cp to all points with coordinates X G f2e,o- This is done in 
the code using the class MappingQEulerian. The material and spatial gradients of an arbitrary field 
can then be interpolated from the material and spatial views of the cell as follows: 

Grad(») ^ ^[Grad^^(X)][(•)^] and grad(») ^ ^[grad$^(x)][(•)^] , 

I I 

where a shape function in the spatial configuration associated with an arbitrary node / is denoted 
By choosing the arbitrary field in the material view as the spatial placement x one obtains JP, 
while choosing the arbitrary field in the spatial view as the material placement X one obtains /. 

The volume Jacobian determinant J and its inverse j are determined from the ratio of the 
Jacobian determinant of the map from the isoparametric configuration to the reference configuration 
J\j,o lo the inverse of the Jacobian determinant of the map from the isoparametric configuration 
(reference cell) to the spatial configuration J\j. 

An identical approach to the volume is followed on the surface. The functionality provided 
by the codimension routines in deal. II allows the surface gradient of a field to be evaluated in a 
straightforward manner. The surface deformation gradient and its inverse are thus computed without 
needing to invert F. This is critical as F is rank deficient. The surface Jacobian determinant 
J and its inverse j are determined from the ratio of the Jacobian determinant of the map from 
the isoparametric configuration to the reference configuration q to the inverse of the Jacobian 
determinant of the map from the isoparametric configuration to the spatial configuration J\j. 
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volume 




Figure 3: The kinematics of the motion in the volume and on the surface. Also shown in the 
Lagrangian and Eulerian view of a cell relative to the reference cell. 

4.4 Parallelization of key tasks 

The library TBB (TBB, 2013) is used to perform as many computationally intensive distributed 
tasks as possible. These include: 

1. the assembly of the tangent and residual contributions from the cells in the volume and on the 
surface; 

2. updating the state of the ContinuumPoint after solving the linearized problem in the Newton 
scheme. 

The main tool for using TBB in deal. II is the WorkStream class. The bottleneck, however, for 
large computations is the linear solver. The code could, in theory, be parallelized in a distributed 
sense relatively easily. The primary challenge would be to ensure that the independent distributed 
volume and surface meshes can communicate, via a map of the common degrees of freedom. 


10 















4.5 The choice of linear solver and preconditioner 

The choice of linear solver and preconditioner are specified in the parameter file parameter .prm. 
The default choice is the conjugate gradient solver (the matrix problem is symmetric) with Jacobi 
preconditioning as provided by the deal. II library. Both this choice of preconditioner and solver 
are multithreaded. The choice works well for the majority of problems investigated. 

5 Numerical results 

The objective of the current section is to elucidate key features of the formulation using four example 
problems. The first example is the Cook’s membrane problem. This is a widely-used benchmark. 
The purpose of the example is to provide detailed information for a series of uniformly refined meshes. 
Neo-Hookean type surface effects are considered. The second example illustrates neo-Hookean type 
surface effects in a nanowire undergoing significant tensile extension. Surface tension is omitted, 
i.e. 7 = 0. The third example of a liquid bridge explores the role of isotropic surface tension in a 
membrane surrounding a thin-walled cylinder with vanishing material properties. This example is 
a good benchmark for membrane problems as an analytical solution is available. It is however not 
ideally suited to problems in surface elasticity as, by definition, we require the energetic surface to 
be the boundary of a volume. The fourth example illustrates neo-Hookean type surface effects in 
a nanoscale plate with a realistic rough surface. This example is used to assess the performance of 
the implementation for relatively complex and realistic geometries. 

The constitutive relations are given in Summary 2. Trilinear and bilinear elements are used 
in the volume and on the surface, respectively. This can be modified in the input parameter file. 
These results have been discussed, and additional context provided, in the companion paper (Javili 
et ah, 2013b). The objective here is to provide information on the performance and features of the 
numerical scheme. 

5.1 Cook’s membrane: neo-Hookean boundary potential 

The Cook’s membrane is a widely-used benchmark problem in solid mechanics. Consider the can¬ 
tilever beam shown in Fig. 4. The left face is fully fixed and a uniform traction applied to the right 
face. All of the remaining faces possess a surface energy (i.e. the back and front, and top and bottom 
faces). 

The coarsest triangulation of the volume is obtained by uniformly dividing the horizontal and 
vertical edges into ten segments and creating a (non-afhne) triangulation consisting of 10x10x1 cells. 
Subsequent triangulations are obtained by uniformly refining this initial one. The Neumann traction 
is applied in ten equal step. 

The material properties in the volume are given in Fig. 4. The material properties of the energetic 
surface are varied while maintaining the ratio A//i = A//i = 1.5. The following two norms in the 
volume and on the surface 

[F : P]^dV 

\Bo 

and the magnitude of the displacement of the point midway between the front and back faces on 
the edge common to the right and top faces, labelled A, for various levels of mesh refinement is also 
shown in Fig. 4. 




and 

[[F : P]^ dA 


J 
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o 

II 

II 

Refinement level 
(number of cells) 


\ua\ 



I'^aI 

0 

(10x10x1) 

1.04019e+05 

13.6392 

5.91255e+04 

5.97436e+04 

8.54514 

1 

(20x20x2) 

9.91750e+04 

14.3084 

5.57675e+04 

5.95120e+04 

8.99501 

2 

(40x40x4) 

9.71923e+04 

14.5326 

5.45158e+04 

5.92104e+04 

9.13706 

3 

(80x80x8) 

9.64376e+04 

14.6100 

5.40838e+04 

5.90174e+04 

9.18937 

4 

(160x160x16) 

9.61404e+04 

14.6382 

- 

- 

- 


Figure 4: The initial 10x10x1 triangulation of the material configuration for the Cook’s membrane. 
The fixed material properties in the volume are given. The final spatial configuration for the ratios 
A/A = 0 and A/A = 1 are shown. The measured norms in the volume and on the energetic surface, 
and the displacement of point A, at four levels are refinement are recorded. Note that for A/A = 1 
and a refinement level of 4, the simulation aborted. 


The presence of the energetic surface significantly reduces the amount of deformation. A sub¬ 
stantial amount of energy goes into deforming the energetic surface that would otherwise go into 
deforming the volume. Numerical problems are encountered in the simulation on the finest mesh 
(160x160x16) when the energetic surface is present. The simulation aborts in the final time step 
due to a negative Jacobian determinant at a quadrature point. The reason for this appears to be 
due to the inherent length scale that the finite element implementation of surface elasticity theory 
introduces. The surface area to volume ratio of a cell (assumed to be a cube) scales as 4/h where 
h is the edge length. If all the surfaces of the cell are assumed energetic, the contribution from the 
energetic surfaces of a volume cell will 64 times greater (relative to the volume contributions) for the 
finest mesh relative to the coarsest. This can clearly introduce numerical problems and is similar 
to the issues that arise in heterogeneous material with vastly different material properties. This 
apparent deficiency is currently under investigation. In related work, the influence of the spatial 
discretization in conjunction with the use of negative surface parameters, which have been reported 
in the literature, is assessed in Javili et al. (2012). 
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Although the results are clearly converging upon mesh refinement the presence of an extremely 
high stress concentration in the cells along the upper left edge of the domain limits the convergence 
rate. 

5.2 Nanowire: neo-Hookean boundary potential 

One application of surface elasticity theory is to describe surface effects in nanowires (see e.g. He and 
Lilley, 2008; Yun and Park, 2009; Yvonnet et al., 2011). Consider the benchmark example shown 
in Fig. 5. The front and back pentagonal faces of the wire are prevented from displacing in the X 
and Y directions. The wire is extended in the Z-direction by an amount 2 (i.e. 40% of the original 
length). The unconstrained surfaces on the side of the wire possess a surface energy. 



Figure 5: The triangulation of the material configuration for the nanowire. The fixed material 
properties are given. The final deformed (spatial) configuration of the volume and the surface for 
three different ratios of A/A are shown. 

The triangulation of the volume, provided in nanowire_f ine. inp, is more refined at the inter¬ 
section of the faces on the surface and towards the front and back faces. This is done to better 
resolve the expected stress concentrations. The volume and surface are discretized into 45 570 and 
4500 elements, respectively. The prescribed deformation is applied uniformly in 10 steps. An alter- 
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native strategy to generate a refined mesh would be adaptive mesh refinement. A coarse mesh is 
also provided in the file nanowire_coarse. inp. 

The Lame parameters in the volume are fixed at the values given in Fig. 5. Similarly, the neo- 
Hookean energetic surface is characterised the surface Lame parameters A and /I. The surface 
shear modulus /I is set to zero and the ratio A/A varied. 

The response in the absence of a surface energy is given by choosing A/A = 0 and is shown 
in Fig. 5. The stress in the volume concentrates at the corners on the front and back faces. The 
initially pentagonal cross section reduces uniformly in size along the length. 

A surface energy is then assigned to the external surface. The stress in the volume P concentrates 
along the lines forming the intersections of the external surfaces as the value of A increases. Increasing 
the surface energy causes the resulting deformed cross section to tend to circular, thus increasing 
the stress in the volume in regions where the surface curvature is not smooth. 

The convergence history of the Newton scheme during the first timestep for A/A = 1 is given in 
Listing 1. The various columns under the heading SOLVER STEP correspond to, for each iteration, 
the application of the constraints (CST), assembling the system matrix contributions from the vol¬ 
ume and the surface (ASS_v and ASS_s, respectively), the solution of the linear system (SLV), and 
updating the data stored at the continuum points in the volume and on the surface (UCP_v and 
UCP_s, respectively). The two columns labelled LIN_IT and LIN_RES record the number of itera¬ 
tions required by the linear solver and the final converged residual value, respectively. The final two 

□ 

columns give the norms of the global residual vector R and the increment in the solution Au, defined 
in (5). The computation of these norms is restricted to unconstrained degrees of freedom. The 
Newton procedure continues until these two measures of the error decrease below the user-specified 
tolerance. 


Listing 1: The convergence history of the Newton scheme during the first timestep for the nanowire 
problem. 

Timestep 1 @ 0.1s 




SOLVER STEP 


1 

LIN_IT 

LIN_RES 

1 1 RnORM 1 

1 dUTlORM 1 

0 

CST 

ASS_v 

ASS_s 

SLV 

UCP_v 

UCP_s 1 

274 

1.064e+02 

l.OOOe+00 

l.OOOe+00 

1 

CST 

ASS_v 

ASS_s 

SLV 

UCP_v 

UCP_s 

795 

2.873e+00 

2.715e-02 

3.603e-03 

2 

CST 

ASS_v 

ASS_s 

SLV 

UCP_v 

UCP_s 

864 

4.630e-02 

4.425e-04 

4.877e-05 

3 

CST 

ASS_v 

ASS_s 

SLV 

UCP_v 

UCP_s 

1026 

6.109e-04 

6.134e-06 

2.349e-07 

4 

CST 

ASS_v 

ASS_s 

SLV 

UCP_v 

UCP_s 

1440 

3.414e-08 

3.187e-10 

1.660e-10 


CONVERGED! 


Relative errors: 

Solution: | dUTlORM | 1.660e-10 

Residual: |RnORM| 3.187e-10 


5.3 Liquid bridge: isotropic surface tension effects 

A liquid bridge is a thin film suspended between two rigid side walls, as depicted in Fig. 6. Surface 
tension acts in the material (initial) configuration to deform the surface so as to minimise its surface 
area. In order to model the liquid bridge using the approach to surface elasticity adopted here, a 
volume must be present. To that end, a neo-Hookean thin-walled cylinder (wall thickness = 0.1) 
with Lame parameters of A = 0 and /i = 10 is enclosed by the energetic surface. Ideally, the volume 
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Figure 6 : The liquid bridge in the material configuration (A) (volume and surface) and the final 
deformed state of the surface (B). The deflection of the midpoint of the surface versus the value of 
7 is shown. 


would have vanishing strength but numerically this results in a poorly conditioned system matrix. 
The volume and surface triangulation contains 36 226 and 18 290 elements, respectively. 

The surface free energy, i.e. 7 , is linearly increased over 20 equal steps to a value of 100 and the 
displacement of a point on the middle of the surface monitored. The deflection of a midpoint on 
the surface for varying 7 is shown in Fig. 6 . The majority of the deformation occurs for 7 < 10. 
Thereafter, as 7 increases so the midpoint deflections converges to the analytical solution of 0.6373 
(see Javili and Steinmann, 2010). 

The numerical solution of this problem is somewhat sensitive. The rate at which the solution 
converges to the exact solution with increasing surface tension 7 is governed by the material prop¬ 
erties of the volume, i.e. /i, and the spatial discretization. The correct way to solve this problem 
robustly would be to disregard the volume and consider the membrane alone. This could be done 
relatively easily but would require the structure of the code to be modified significantly. 
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5.4 Bending of a nanoscale plate with a rough surface 

Surface roughness can have a significant, and often complex, influence on the response of nanoscale 
objects. Consider the cantilever plate shown in Fig. 7. The profile of the upper surface was produced 
using an open-source rough surface generation tool (Bergstrom, 2013).^ The rough upper surface 
is assumed energetic. All other surfaces are planar and standard. The Lame parameters in the 
bulk ^e fixed at the values given in Fig. 7. The surface shear modulus is set to zero and the 
ratio A/A varied. The left edge of the plate is fully fixed and a prescribed surface traction of 
6o = is imposed incrementally (in 10 uniform steps) on the lower face. The bulk and 

surface triangulation contains 250 000 and 10 000 elements, respectively (a coarser triangulation is 
provided in the accompanying online documentation). 



Figure 7: The triangulation and material configuration for the cantilever plate. The fixed material 
properties are given. The final spatial configuration of the bulk for two different ratios of A/A are 
shown. The magnitude of the displacement u is plotted. 

The response of the plate to the loading for A/A = 0 and A/A = 1 is compared in Fig. 7. The 

^The routine used, rsgene2D, produces a Gaussian height distribution with an exponential auto-covariance. The 
input parameters were 100 divisions, a surface length of 2, a root mean square height of 0.05, and an (isotropic) 
correlation length of 0.25. The surface was then scaled uniformly to have a surface length of 8. 
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plate with the energetic rough surface deforms considerably less. 

The example of the cantilever plate with the rough surface provides a good test of the robust¬ 
ness and efficiency of the numerical implementation. The Newton scheme exhibited the quadratic 
convergence associated with a consistently derived tangent. A summary of the time spent in the 
various sections of the code is given in Listing 2. It is clear that for problems of a reasonable size, 
the bottleneck is the solution of the linear system which accounts here for 77.7% of the simulation 
time. A distributed parallel implementation utilising an AMG preconditioner, in the spirit of the 
approach used by Frohne et al. (2013), would be one way to address this issue. 


Listing 2: Summary of the time spent in the various parts of the code for the cantilever plate 
problem. The simulation was performed on a dual processor 8-core Intel Xeon 2.4GHz GPU with 
64GB RAM 


1 

1 

Total wallclock time elapsed 

since start 

1 

1 

6.375e+03s 

+ 

1 

1 


-+ 

1 

1 

1 

1 

1 

Section 

1 no . calls 

1 

1 

1 

1 

wall time 

1 

1 

1 

% of total 

1 

1 

1 

1 

1 

Assemble system volume 

1 

1 55 

1 

1 

2.667e+02s 

1 

1 

4.18e+00% 

1 

1 

1 

Assemble tangent surface 

55 

1 

1.795e+00s 

1 

2.8 2e-02% 

1 

1 

Construct grid 

1 

1 

7.110e+00s 

1 

1.12e-01% 

1 

1 

Linear solver 

55 

1 

4.950 e+03s 

1 

7.76e+01% 

1 

1 

Postprocess results 

11 

1 

8.970eT02s 

1 

1.41 e+01% 

1 

1 

Setup system 

1 

1 

8.555e+00s 

1 

1.34e-01% 

1 

1 

Update CP data 

110 

-H- 

1 

2.083e+02s 

1 

+ 

3.2 7e+00% 

1 

-+ 
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Summary 2 Hyperelastic material models for the volume and surface 

Volume: 


J := DetJP : Jacobian determinant 


/i , A : Lame constants 


Free energy 


^(F) = I Aln^ J + i/i[F : F- 3 - 21 n J] 


Piola-Kirchhoff stress 


P{F) = — = \hiJf + l,[F-f] 


Piola stress tangent 


BP 

HF) = ^ = A [/‘ 0 /‘ + In J D] + M [I - D] 


dF 

D:= — = - 

OF ■' 


^ dF 


Surface: 

J := DetF : Surface Jacobian determinant yU, A : Surface Lame constants 7 : Surface tension 


Surface Free energy 


^'(F) = iAln^ J + i/x[F : F- 2 - 21 nJ] + 7 J 


Surface Piola-Kirchhoff stress P{F) = = A In J + /I [F — /^] + 7 J 

BF 

^ ^ Qp __ __ 

Surface Piola stress tangent A{F) = = A [/^ 0 + In J D] + /I [I — D] + 7 J [/^ (g) + D] 




7 dF .-7 
I := = ^(g) J 

BF 
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6 Conclusion 


An efficient finite element scheme for the solution of problems in nonlinear surface elasticity was 
presented. The effectiveness of the scheme was demonstrated using three example problems. The 
complete, documented source code has been provided. The code also provides a template for prob¬ 
lems in classical finite elasticity. 

Various extensions of the approach are planned. The first is to perform adaptive mesh refinement 
using the tools provided by deal .II. The key challenge here is the projection of the data conveniently 
located at the level of the quadrature point between successive meshes. This should not be too 
problematic as, unlike plasticity, none of the governing (evolution) equations are defined pointwise. 

The second extension would be to improve the efficiency of the code by implementing distributed 
parallelism. The primary challenge would be to ensure that the independent distributed volume 
and surface meshes can communicate, via a map of the common degrees of freedom. The use of 
better-suited parallel solvers and preconditioners is also under investigation. Based on the results 
presented in Frohne et al. (2013), for similar problems in elastoplasticity, the choice of the AMG 
preconditioner and BiCGStab solver should provide a good starting point. 

The tangent is recomputed during each iteration of the Newton scheme. This operation is not 
that expensive, relative to the solution of the linear system, but may not be necessary. Finally, the 
Newton scheme should be extended to include a line search algorithm and adaptive time stepping 
introduced. 

A way to handle the numerical problems that arise as the (energetic) surface area to volume 
ratio increases upon mesh refinement is currently under investigation. 
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Appendix A Library information and source code 

The complete, documented source code can be found online at www.cerecam.uct.ac.za/code/ 
surf ace_energy/doc/htmI/index .html. The code was compiled with version 8.0.0 of the deal. II 
library using the precompiled Mac OS X package, and on a Linux machine using version S.l.pre. 

The steps to follow to run the numerical examples are described in the online documentation. 
The procedure to compile the code is identical to the deal. II tutorials. 


Appendix B Features of the implementation in deal. II 

The structure of the implementation follows the majority of the deal. II tutorial examples. The 
following section explains the key namespaces, classes and methods in the implementation. 

B.l Namespaces 
namespace Surface_EIasticity 

The complete implementation is wrapped within this namespace. 
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namespace Parameters 

The parsing and reading of the input parameters from the various parameter files are handled in 
this namespace. 

namespace AdditionalTools 

This namespace defines several operations from tensor algebra that are not part of the deal. II 
library. 

B.2 Classes and methods 
class Time 

Simple time management class used to incrementally advance the quasi-static problem. 

class ContinuumPoint 

As discussed in Section 4.2, this class is responsible for the handling the constitutive response. This 
class would need to be modified if a different type of material model was adopted. 

class Solid 

This class contains the core routines to manage to the code. A Solid encapsulates the volume and 
the surface and all the routines required to control the problem. The key methods within the class 
are: 

run 0 

This is the method responsible for the high level control. It calls the routines to create the grid, 
solve the problem using the Newton scheme, and to output the results. 

system_setup () 

This method is used to initialise the system matrix and various other vectors, including the incre¬ 
mental solution and the right-hand side. 

assemble_system_volume () 

This method is responsible for assembling the contributions from the volume to the system matrix 
and the right-hand side vector. 

assemble_system_surface () 

This method is responsible for assembling the contributions from the surface to the system matrix 
and the right-hand side vector. 

make_constraints (...) 

The Dirichlet constraints are assembled in this routine. 
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solve_nonlinear_timestep (...) 

This routine controls the Newton-Raphson procedure used to solve the problem. It continues the 
iterative procedure until a converged solution is obtained. 

solve_linear_system 

This routines solves the linear system of equations using the choice of solver and preconditioner 
specified in the parameter file. 

output .results () 

This routine is responsible for outputting the results in vtu format for subsequent post-processing. 
The user can choose to write the data associated with the quadrature points. Note, this can be 
time consuming for large problems. The solution in the volume and on the surface are outputted in 
separate files. 

update _volume_cp_incremental 

Given an updated solution, this routine is responsible for updating the data stored at the quadrature 
points in the volume. 

update.surf ace_cp_incremental 

Given an updated solution, this routine is responsible for updating the data stored at the quadrature 
points on the surface. 


References 

2013. URL www.threadingbuildingblocks.org. 

W. Bangerth, R. Hartmann, and G. Kanschat. deal.II — a general purpose object oriented finite 
element library. ACM Transactions on Mathematical Software, 33(4):24, 2007. 

W. Bangerth, T. Heister, L. Heltai, G. Kanschat, M. Kronbichler, M. Maier, B. Turcksin, and T.D. 
Young. The deal.II Library, Version 8.0. arXiv.org, page arXiv: 1312.2266, 2013a. 

W. Bangerth, T. Heister, G. Kanschat, et al. deal. II Differential Equations Analysis Library, 
Technical Reference, 2013b. URL http://www.dealii.org. 

D. Bergstrom, 2013. URLwww.mysimlabs.com/surface_generation.html. 

A. DeSimone, L. Heltai, and G. Manigrasso. Tools for the solution of pdes defined on 
curved manifolds with deal.ii. 2009. URL www.dealii .org/7.3.0/reports/codimension-one/ 
desimone-heltai-manigrasso.pdf. 

H.L. Duan, J. Wang, and B.L. Karihaloo. Theory of elasticity at the nonoscale. Advances in applied 
mechanics, 42:1-68, 2009. 

J. Frohne, T. Heister, and W. Bangerth. Efficient numerical methods for the large-scale, parallel 
solution of elastoplastic contact problems. 2013. In review. 


21 



M.E. Gurtin and A. Murdoch. A continuum theory of elastic material surfaces. Archive for Rational 
Mechanics and Analysis, 57(4):291-323, 1975. 

J. He and C.M. Lilley. Surface Effect on the Elastic Behavior of Static Bending Nanowires. Nano 
Letters, 8(7):1798-1802, 2008. 

L. Heltai. On the stability of the finite element immersed boundary method. Computers & Structures, 
86(7-8):598-617, 2008. 

A. Javili and P. Steinmann. A finite element framework for continua with boundary energies. Part 
II: The three-dimensional case. Computer Methods in Applied Mechanics and Engineering, 199 
(9-12):755-765, 2010. 

A. Javili, A. McBride, P. Steinmann, and B.D. Reddy. Relationships between the admissible range 
of surface material parameters and stability of linearly elastic bodies. Philosophical Magazine, 92: 
3540-3563, 2012. 

A. Javili, A. McBride, and P. Steinmann. Thermo mechanics of Solids with Lower-Dimensional 
Energetics: On the Importance of Surface, Interface and Curve Structures at the Nanoscale. A 
Unifying Review. Applied Mechanics Reviews, 65(1):010802, 2013a. 

A. Javili, A. McBride, P. Steinmann, and B.D. Reddy. A curvilinear-coordinate-based finite element 
methodology. A unified computational framework for bulk and surface elasticity theory. In review, 
2013b. 

J. Weissmuller, H.-L. Duan, and D. Earkas. Deformation of solids with nanoscale pores by the action 
of capillary forces. Acta Materialia, 58(1): 1-13, 2010. ISSN 13596454. 

G. Yun and H.S. Park. Surface stress effects on the bending properties of fee metal nanowires. 
Physical Review B, 79(19):32-35, 2009. 

J. Yvonnet, A. Mitrushchenkov, G. Chambaud, and Q.-C. He. Einite element model of ionic 
nanowires with size-dependent mechanical properties determined by ab initio calculations. Com¬ 
puter Methods in Applied Mechanics and Engineering, 200(5-8):614-625, 2011. 


22 



